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We present results from quantum Monte Carlo simulations of trapped bosons in optical 
lattices, focusing on the crossover from a gas of softcore bosons to a Tonks-Girardeau gas in a 
one- dimensional optical lattice. We find that depending on the quantity being measured, the 
behavior found in the Tonks-Girardeau regime is observed already at relatively small values of 
the interaction strength. A finite critical value for entering the Tonks-Girardeau regime does 
not exist. Furthermore, we discuss the computational efficiency of two quantum Monte Carlo 
methods to simulate large scale trapped bosonic systems: directed loops in stochastic series 
expansions and the worm algorithm. 
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1. Introduction 

Interest in the properties, phase transitions and vari- 
ous phases of strongly correlated systems in reduced di- 
mensionality has a very long history in condensed matter 
physics. The experimental realization of Bose-Einstein 
condensates (BEC) in traps 1 offered the possibility of 
studying such phenomena in precisely engineered sys- 
tems. This goal was eventually achieved by placing the 
condensates on optical lattices 2,3 which has the effect 
of packing the bosonic atoms in close proximity to one 
another and even producing multiple occupancy of lat- 
tice sites. Consequently, strongly correlated systems are 
produced on such optical lattices which have the advan- 
tage of being defect free. However, the lattices are im- 
mersed in the confining trap and are consequently not 
translationally invariant. This makes delicate the task of 
applying results known in condensed matter physics for 
translationally invariant systems to BEC on optical lat- 
tices. The rush to do this has led to some misunderstand- 
ings especially concerning "quantum phase transitions" 
on the optical lattice. 

Our goal here is to use quantum Monte Carlo (QMC) 
methods, such as the stochastic series expansion (SSE) 
and the worm algorithm, to study in some detail the 
properties of BEC on optical lattices and elucidate the 
nature of the observed "transition" . 

The paper is organized as follows. In section 2, we 
introduce the bosonic Hubbard model used for the de- 
scription of a cold atomic Bose gas trapped on optical 
lattices and discuss and compare the techniques to sim- 
ulate it. We then review the results presented in Ref. 4 
in section 3, focusing on observables particularly suited 
for the detection of local "phases" or structures and the 
absence of quantum criticality in these inhomogeneous 
systems. In section 4, new results on the transition to a 
Tonks-Girardeau gas regime for bosonic systems in one- 
dimension are presented. 



2. Quantum Monte Carlo simulations 

2.1 Model 

We consider the following Bose Hubbard Hamiltonian 
-tJ2(4bj + h.c.) (1) 

(i,3) 
U 

+ 2 



H 



1) 



as a good description of the low energy physics of cold 
confined bosonic gases. 5 Here t is the nearest neighbor 
hopping, U the onsitc repulsion between the softcore 
bosons, fi the chemical potential and V the curvature 
of the parabolic confining potential imposed by the trap. 
We use standard notations for the bosonic operators. We 
furthermore define /x° ff = fi— Vrf, which is the local effec- 
tive chemical potential experienced by a boson at site i. 

2.2 Quantum Monte Carlo 

The Hamiltonian (1) is particularly well suited for sim- 
ulations using QMC techniques. Indeed and to our best 
knowledge, such techniques are the only ones able to sim- 
ulate efficiently the Hamiltonian (1) for the following rea- 
sons : 

• They work in any dimension (unlike the density 
matrix renormalization group (DMRG) which is re- 
stricted to one-dimensional systems). 

• They can be used to obtain both finite and zero 
temperature properties. 

• They allow the simulation of non homogeneous trap- 
ping potentials (unlike series expansion). 

• Large system sizes can be reached (unlike using ex- 
act diagonalization) . 

• They perform an exact treatment of all terms in the 
Hamiltonian (1) (unlike most analytic approaches). 
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2.3 Efficiency of non-local quantum Monte Carlo 
methods 

We used the most recent world-line QMC algorithms, 
the SSE algorithm 6 with directed loops 7,8 and the worm 
algorithm 9 to perform the simulations based on the open 
source implementations of the ALPS project. 10 

The SSE representation 6 starts from a Taylor expan- 
sion of the partition function in orders of the inverse 
temperature (3: 

00 on 

Z = Trexp(-^) = V ^Tr(-ff)" (2) 

71=0 

00 on 

= E^ E {h\-H\iz)---{i n \-H\h). 

n=0 ' {i 1 ,...i n } 

This representation of the partition function does not 
suffer from errors introduced by time discretization and 
the sampling algorithm can be combined with non-local 
updates such as directed loops. 7 ' 8 However, the SSE 
representation of the partition function corresponds to 
a perturbation expansion in both, diagonal and off- 
diagonal, terms of the Hamiltonian. For the simulation 
of trapped bosonic systems we are naturally interested 
in the limit where the diagonal matrix elements (e.g. the 
trapping potential) are large. In this limit the SSE sam- 
pling can face a considerable slowdown as the probabil- 
ity of exchanging off-diagonal and diagonal bond terms 
in the Hamiltonian is suppressed. On the contrary, the 
worm algorithm 9 samples world lines in the path integral 
representation of the partition function which treats only 
the off-diagonal terms in the Hamiltonian as a perturba- 
tion 

Z = Tr(e-^) =Tr(e-^Te-^ v M) . 

(3) 

The directed loop and worm algorithms both work in 
an extended configuration space, which in addition to 
closed world line configurations, allows for the presence 
of an open world line fragment, the "worm" or "directed 
loop", which is formally introduced by adding a source 
term to the Hamiltonian 

H m = H-r ) J2(bi+bi) ■ (4) 

This source term allows the world lines to be broken 
with a matrix element proportional to n. To generate 
new closed loop configurations in the sampling process, a 
worm is created and a random sequence of local updates 
of the worm is performed where each move fulfills de- 
tailed balance. Nevertheless, this update procedure can 
perform non-local changes of the world line configura- 
tions as the worm can wind around the lattice in the 
temporal or spatial direction and thereby change the par- 
ticle or winding number respectively. 

To discuss the computational performance of the 
two QMC algorithms in simulating confined ultra-cold 
bosonic atoms in two dimensions we compare the sta- 
tistical error obtained for the Monte Carlo estimates of 
observables from simulation runs with a fixed amount 
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Fig. 1. Comparison of the statistical error Aft local of the local 
compressibility re local of a 2D trap with 35 X 35 sites obtained 
from simulation runs with fixed CPU-time for stochastic series 
expansion (SSE) quantum Monte Carlo and the worm algorithm. 
The error bars are obtained by averaging over estimates with 
identical radial distance. The inset shows the ratio of statistical 
errors Ak (SSE) /AK local (worm) for the two algorithms using 
the open-source implementations of the ALPS project. 10 



of CPU time (typically around 100 CPU hours on an 
Opteron 1.8 GHz processor). 

While in the homogenous Bose-Hubbard Hamiltonian 
without confining potential the SSE representation is 
found to perform better, except from the extreme soft- 
core case, 11 the situation is different when simulating 
harmonic traps. In Fig. 1 results of our simulations are 
shown for a two-dimensional system with 35 x 35 sites 
and trap parameters which allow for a Mott plateau 
with integer density ((n^) = 1) in the center of the trap 
(R/L < 0.2). Before discussing the physics of this system 
in the next sections we first consider the statistical er- 
ror of the local compressibility K local (defined by Eq. (5) , 
below), which is shown in Fig. 1 as a function of the dis- 
tance from the trap center. While for both algorithms 
the statistical error is largest for the supcrfluid ring sur- 
rounding the central Mott plateau (0.2 < R < 0.4), the 
worm algorithm yields a statistical error which is a factor 
of 3 smaller than the one obtained for the SSE estimate 
with the same amount of CPU-time. The superior per- 
formance of the worm algorithm is a consequence of the 
small ratio of off-diagonal to diagonal matrix elements 
for the used value of t/U — 0.03 in the simulated Bose- 
Hubbard Hamiltonian. Further reducing this ratio, e.g. to 
study trapped bosonic systems which exhibit more than 
one Mott plateau, will result in an even larger speedup 
of the worm algorithm. 

Fig. 2 shows the ratio of statistical errors of QMC 
estimates for three different observables obtained using 
the SSE and worm algorithm. The statistical errors are 
shown versus the number of lattice sites. For the particle 
number n and the compressibility K local (c.f. Eq. (5)) spa- 
tial averages are shown. There appears to be a constant 
improvement of the worm code over the SSE algorithm, 
independent of the system size, implying that both algo- 
rithms exhibit the same scaling in system size. From our 
simulations we conclude that the worm algorithm is the 
superior QMC technique to simulate trapped ultra-cold 
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Fig. 2. Scaling of the speedup of the worm algorithm compared 
to the SSE algorithm versus system size for the same trap pa- 
rameters as in Fig. 1. Shown are results for the spatially averaged 
local density n and compressibility K local as well as the energy 
E. For all simulations the open-source implementations of the 
ALPS project 10 were used. 




Fig. 3. Spatial dependence of the local compressibility k 1oc , 
of bosons in a two dimensional parabolic trap with curvature 
V/U = 0.002, for fi/U = 0.37 and U/t = 25. A superfluid ring 
surrounding a n = 1 central Mott plateau is clearly resolved. 



bosonic systems. 

3. Simulations of trapped bosonic atoms 

3.1 Coexistence of phases and identification of local 
structures 

As shown in various studies (see for example Refs. 2,3 
for experimental and Refs. 5, 12 for theoretical/numerical 
investigations), the inhomogeneity in the system due to 
the trapping potential induces a co-existence of super- 
fluid and Mott insulator- like regions. This can be un- 
derstood on qualitative grounds by looking at density 
profiles: a region of space with constant integer density 
of particles is interpreted as a Mott insulating region 
(a Mott plateau) , while non-integer densities correspond 
to superfluid regions. To distinguish the state near a 
given site more precisely, we proposed a more quanti- 
tative probe, which will be reviewed here. 4 

For the homogeneous case, the Mott insulating phase 
has a vanishing compressibility k. 13 One way of locally 
characterizing a region of space in the inhomogeneous 
system is to study the local compressibility at a given 
lattice-site i: 

^° Cal = (c^) =/3((n2n) ~ (n ' i)(n)) ' (5) 

which expresses the total density response of the sys- 
tem to a local change of the chemical potential at site i. 
This local compressibility directly corresponds to the to- 
tal compressibility n in the homogeneous case, and serves 
to distinguish the local states in the inhomogeneous sys- 
tem: K l ° cal is zero in a Mott insulating region while it 
remains finite superfluid region. In Fig. 3, we show the 
spatial distribution of K' ocal in a two dimensional trapped 
system. Regions of space with vanishing «;' ocal are clearly 
identified and correspond to Mott insulating regions in 
the sample. 

For the Hamiltonian in Eq. (1), the kinetic term t, 
the repulsion U and the trap curvature V are tunable 
parameters in the experiments, 2 ' 3, 19-22 and a "quantum 




U/t 

Fig. 4. Evolution of the coherence fraction (the height of the 
peak, n(0, 0)) and of the full width at half maximum (FWHM) 
of the coherence peak as a function of U /t for bosons in a two di- 
mensional parabolic trap with curvature V/U = 0.002 and with 
fi/U = 0.37. The threshold for Mott plateau formation is indi- 
cated by the dashed line. 



phase transition" between a Mott insulator and a super- 
fluid can be observed by changing e.g. the value of U/t. 
To study this transition, we calculated the momentum 
distribution function 

"(k) = ^E e ' (r *" r3 ' )k ( 6 l^)> ( 6 ) 

where N is the total number of particles in the system. 
The momentum distribution in Eq. (6) is normalized, and 
the coherence fraction is given by the height of the coher- 
ence peak, n(k — 0). Not only is this quantity useful to 
detect the formation of a Mott plateau as will be shown 
below: it also is a quantity of interest for experiments, as 
n(k) can be determined from interference patterns. 14 

Fig. 4 displays both the coherence fraction n(k = 0) 
and the full width at half maximum (FWHM) of the 
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Fig. 5. Local density, n, and local compressibility, k}° , for a two 
dimensional parabolic trap on a square lattice and for a ladder, 
as functions of fi cS . The parameters for the trapping potential 
are V/U = 0.002, at fi/U = 0.37, and the parameters of the 
ladder model are chosen as to cover the whole superfluid region 
(see Ref. 4). 



Fig. 6. Local density, n, and local compressibility, K local , for the 
Bose-Hubbard model on ladders with different lengths as a func- 
tion of /i eff . The other parameters of the ladder model are chosen 
as in Fig. 5. The inset shows the compressibility for the uniform 
onc-dimensional case as a function of fi for a chain with 64 sites, 
for the same value of U/t as used in the ladder model. 



coherence peak in a two dimensional trap as obtained 
from the QMC simulations. The ratio U/t at which a 
Mott plateau forms in the center of the trap is marked 
by a vertical line. Whereas the coherence fraction does 
not show any specific feature at the threshold for Mott 
plateau formation, we observe a change of the curvature 
in the FWHM at this point. It was shown, 4 that this ap- 
pears as a generic feature and that the FWHM can thus 
serve as an indicator of Mott plateau formation in a con- 
fined Bose gas. We also showed that the appearance of 
fine structure in the momentum distribution (secondary 
peaks) is not related to the Mott plateau formation, dis- 
proving claims in previous work. 14 

3.2 Absence of quantum criticality 

Fig. 3 suggests that the superfluid region identified 
by a non vanishing local compressibility is confined to 
a shell around the central Mott-insulating region. One 
might then expect that in a two-dimensional trap this 
shell essentially behaves like a one-dimensional bosonic 
chain, and displays similar quantum critical behavior. 
For example, the compressibility of a one-dimensional 
bosonic chain diverges at the superfluid-insulator transi- 
tion. 15 However, simulations in a trap indicate no critical 
features in quantities like the local compressibility. 4 

To further investigate this phenomenon and to sim- 
plify simulations, an effective bosonic ladder model was 
introduced in Ref. 4, that represents the superfluid ring 
in the original two-dimensional trap. In this model, each 
leg of the ladder is assigned a different chemical poten- 
tial, which is chosen such that the first leg is always in 
a (n) = Mott insulating phase and the last leg in the 
(n) = 1 Mott insulating phase. The chemical potential 
is then interpolated (for example linearly) between these 
two values for the other legs of the ladder. 

In Fig. 5, we show the local density and the local com- 
pressibility K local as a function of the effective chemical 
potential /j, for both simulations, the 2d trap and the 
ladder model. The overall good agreement of the two sets 



of curves indicates that indeed the ladder model captures 
the essential physics of the trapped two-dimensional sys- 
tem. The density curves from the trap and the ladder 
model coincide almost perfectly, indicating the validity 
of a local density approximation. 4 There are small differ- 
ences in the local compressibility between the trap and 
the ladder system, which are due to the different shape 
of the local chemical potential in the two systems. Nev- 
ertheless, the two different compressibility curves share 
the same overall shape. 

Having shown that the ladder model describes cor- 
rectly the different physical features found in the trapped 
system, we can study the occurrence of quantum crit- 
icality in this simpler model. In Fig. 6, the local den- 
sity and the local compressibility of the ladder model 
are shown for two different system sizes. We observe two 
very broad peaks in the compressibility, which might in- 
dicate quantum critical behavior. However, these peaks 
do not sharpen upon increasing the length of the ladder, 
as seen in Fig. 6. The inset of Fig. 6 displays the com- 
pressibility of a one-dimensional chain which, in contrast 
with the ladder results, displays very sharp peaks at the 
insulating-superfluid phase transition, even for a system 
of moderate size (L = 64). 

We thus find that no one-dimensional quantum critical 
behavior is found in the realistically simulated trapped 
systems. This is in agreement with the absence of "crit- 
ical slowing down" observed in the experiments. 2 ' 3 We 
thus interpret the "quantum phase transition" seen in the 
experiments as a crossover of changing fractions of the 
Mott-insulating and superfluid regions. Critical behavior 
could however be observed for flat confining potentials. 
For a more elaborate discussion on the ladder model and 
the absence of quantum criticality, we refer the reader to 
Ref. 4. 
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4. Approaching the Tonks-Girardeau regime 

4-1 The Tonks-Girardeau gas 

In previous sections, we studied properties of con- 
fined ultra-cold bosonic atoms in optical lattices using 
the Bosc-Hubbard Hamiltonian with an on-site repulsion 
U/t. In the limit of infinite repulsion, U/t — > oo, this 
model maps onto a model of hardcore bosons which are 
constrained to a maximum occupancy of one boson per 
site. In the continuum limit of the one-dimensional Bosc 
gas, the regime of infinite repulsion defines the Tonks- 
Girardeau gas, 16, 17 which in many respects behaves sim- 
ilar to fermions due to the enforced impenetrability. The 
realization of a Tonks-Girardeau gas in trapped bosonic 
atom systems 18 initially turned out difficult, due to the 
restricted interaction range accessible from the exper- 
imental setup. 19 ' 20 However, recently Tonks-Girardeau 
gas behavior was observed in one-dimensional Bose gases 
in the presence of an optical lattice potential. 21 The un- 
derlying lattice structure increases the effective mass of 
the bosons, thus increasing the ratio 7 = -Emt/^kin, of 
interaction energy {E- m t) to kinetic energy (-Ekin)- 21 In 
particular, excellent agreement between the experimen- 
tal data and calculations based on the fermionization 
approach to the Tonks-Girardeau gas was obtained in 
Ref. 21, already for values of 7 « 5. In another recent 
study, 22 observation of Tonks-Girardeau gas behavior 
was reported for similar value of 7 « 5 also in the absence 
of an optical lattice. 

4-2 When do we observe Tonks-Girardeau gas behavior? 

In the following, we compare both the ground state 
density distribution and the momentum distribution 
function of bosons in one-dimensional optical lattices 
with a harmonic confinement potential to address the 
question when a system of repulsively interacting bosons 
shows behavior similar to a Tonks-Girardeau gas. We 
performed simulations for chains of N s — 50 lattice sites, 
typical for experimental setup, 21 and considered different 
fillings of the system by adjusting the chemical potential. 

An appropriate means of quantifying the relevance of 
the interactions in the presence of an optical lattice is the 
ratio U/t, given in terms of the microscopic model pa- 
rameters. 23 To a given value of U/t, an associated effec- 
tive value of 7cff = U / (2t(n)) is obtained by determining 
the particle density (n) during the simulation. 24 

In Fig. 7 the QMC results for j e s are shown as a func- 
tion of U/t for \i/t = and 0.3. For a given value oiU/t, 
the effective value of the interactions, 7 e fj, is larger for 
/j,/t = than for fi/t = 0.3, due to a lower bosonic density 
in the first case. We thus expect the behavior of softcore 
bosons to fit well to that of a Tonks-Girardeau gas for 
smaller values of U/t for fi/t — than for fi/t = 0.3. 

In Fig. 8 (a,b) the spatial density distribution , n{x), 
and the momentum distribution function, n(k), are 
shown in the low-density regime for fi/t = 0, with no 
Mott-insulating region present in the trap. We find that 
for U/t = 50 both quantities are already close to the 
hardcore limit. Data taken at U/t = 100 (not shown) did 
not exhibit any visible difference to the hardcore limit in 
this regime. In contrast, for U = 5 the density in the trap 
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Fig. 7. Effective interaction strength 7 e ff as a function of the ra- 
tio U/t for bosons on one-dimensional optical lattices , inside a 
harmonic confinement potential of strength V/t = 10~ 3 . Results 
of quantum Monte Carlo simulations arc shown for different val- 
ues of the chemical potential, fi/t = 0, and 0.3. The inset shows 
the average deviation e n in the local density between the softcore 
and hardcore cases as a function of the interaction strength U/t 
for fi/t = and 0.3. 
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Fig. 8. Density distribution, n(x), and momentum distribution 
function, n(k), in the ground state of softcore bosons on one- 
dimcnsional optical lattices for different values of U/t = 5, 10, 
50, and for hardcore bosons (U/t = 00), inside a harmonic con- 
finement potential of strength V/t = 10~ 3 and for fx/t = 0. 



center is significantly larger than in the hardcore case, as 
is the coherent fraction n(k = 0). One might thus con- 
clude, as in Ref. 25, that large values of 7 > 100 are 
needed before the behavior of softcore bosons fits well to 
that of a Tonks-Girardeau gas. Observing that near the 
boundary of the system, the density distributions is close 
in all cases and that due to the low density of bosons in 
this regime, the value of U is irrelevant, we prefer a dif- 
ferent look at this problem. 

For a more detailed analysis, we plot the averaged 
squared deviation e n (U) = l/N s J2i( n i(U) — rii(U = 
00)) 2 as a function of U/t in the inset of Fig. 7. We find, 
that e n (U) decreases algebraically upon increasing U, i.e. 
e n (U) cx U~ a , with a « 2. This shows that the crossover 
of softcore bosons into the hardcore limit is smooth, and 
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Fig. 9. Density distribution, n(x), and momentum distribution 
function, n(k), in the ground state of softcore bosons on one- 
dimensional optical lattices for different values of U/t = 5, 10, 
50, and for hardcore bosons (U/t = oo), inside a harmonic con- 
finement potential of strength V/t = 10~ 3 and for fi/t = 0.3. 



that no finite critical value for U /t exists. 

Considering the momentum distribution function n{k) 
shown in Fig. 8 (b), we find that they are similar 
to the hardcore limit for all values of U/t in this 
intermediate-density regime. The main differences occur 
for the plateau value near k = 0, whereas the character- 
istic slopes of the curves around k/a « ir/2 21 are insen- 
sitive to the value of U/t. Due to this similar structure, 
and since the overall scale of n(k) is not fixed by the ex- 
perimental data, the measured momentum profiles could 
be fitted in Rcf. 21 to the momentum profile of hardcore 
bosons already for moderate values of 7. Experimentally 
it is thus difficult to distinguish the softcore regime from 
the hardcore limit in the low-density region, and behav- 
ior similar to a Tonks-Girardeau gas emerges already for 
moderate values of 7 > 5, as observed in Rcf. 21. 

Next, we consider the case of an increased bosonic den- 
sity inside the trap, obtained by increasing the chemical 
potential. In particular, we consider (j,/t = 0.3, and show 
results for the local density and momentum distribution 
function in Fig. 9. As seen from Fig. 9 (a), hardcore 
bosons develop an extended Mott plateau in the trap 
center for this value of \i. Plateau formation is observed 
for softcore bosons only for values of U/t > 10. Further- 
more, the algebraic decrease in e n (U) upon increasing U 
found for [i/t = is also observed for [i/t = 0.3 and 
with the same exponent of a s=s 2. Turning to the mo- 
mentum distribution functions, we find from Fig. 9 (b), 
that the momentum profile for U/t = 5 and 10 now show 
significant quantitative differences to the one in the hard- 
core limit. In particular, the characteristic slopes of the 
curves near k/a « tt/2 now differ significantly from the 
hardcore case. For high bosonic densities, the differences 
between softcore and hardcore behavior are thus more 
pronounced, and easier detectable by measuring the ex- 
perimentally accessible momentum profile. 

5. Conclusion 

The main new results of the present paper concern the 
crossover of softcore bosons into the Tonks-Girardeau 



regime. From our simulations we conclude that this 
crossover appears smooth, and that no finite critical 
value of the interaction strength U/t exists. Instead, 
upon increasing U/t, differences in the local density dis- 
tributions between the softcore and hardcore case de- 
crease algebraically, with an exponent that is insensitive 
to the density inside the trap. 

In the low-density region, the shape momentum pro- 
files of softcore bosons are qualitatively very similar to 
those observed in the Tonks-Girardeau gas already for 
moderate values of U/t. This explains the observation in 
Ref. 21, where the measured momentum profiles could 
well be fitted using a fcrmionization approach to the 
Tonks-Girardeau gas for 7 > 5. The value of 7 or U/t be- 
yond which an experiment is well described by a Tonks- 
Girardeau gas picture depends sensitively on the quan- 
tity being measured, the measurement error and the den- 
sity of bosons. In particular in the low-density region, we 
find the momentum profile to be rather insensitive to the 
value of U/t, allowing a description of the softcore bo- 
son system in terms of a Tonks-Girardeau gas already 
at small values of 7. For larger densities, the additional 
confinement of the superfluid in the hardcore limit leads 
to quantitative differences in the momentum profile be- 
tween softcore and hardcore bosons and thus larger val- 
ues of U/t are needed in order to observe the behavior of 
a Tonks-Girardeau gas. 
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